sf can represent all common vector geometry types (raster data classes are not supported by sf): points, lines, polygons and their respective ‘multi’ versions (which group together features of the same type into a single feature). sf also supports geometry collections, which can contain multiple geometry types in a single object.

Source: Geocomputation with R by Robin Lovelace

Coordinate reference system (CRS)

How the spatial elements of the data relate to the surface of the Earth. That coordinate systems are a key component of geographic objects

Knowing which CRS your data is in, and whether it is in geographic (lon/lat) or projected (typically meters), is important and has consequences for how R handles spatial and geometry operations

Source: Geocomputation with R by Robin Lovelace

Load Swiss geo data

Data from Office fédéral de la statistique - Limites communales généralisées: géodonnées

muni <- st_read("shp/g2g22.shp")

## Reading layer `g2g22' from data source 
##   `/Users/duc-q.nguyen/Documents/ddj/HEIGVD/R/shp/g2g22.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2148 features and 19 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 2485424 ymin: 1075268 xmax: 2833837 ymax: 1295937
## Projected CRS: CH1903+ / LV95

Anatomy of a sf object

str(muni)

## Classes 'sf' and 'data.frame':   2148 obs. of  20 variables:
##  $ GMDHISTID: num  13256 11742 11801 11992 12249 ...
##  $ GMDNR    : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ GMDNAME  : chr  "Aeugst am Albis" "Affoltern am Albis" "Bonstetten" "Hausen am Albis" ...
##  $ BZHISTID : num  10053 10053 10053 10053 10053 ...
##  $ BZNR     : num  101 101 101 101 101 101 101 101 101 101 ...
##  $ KTNR     : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ GRNR     : num  4 4 4 4 4 4 4 4 4 4 ...
##  $ AREA_HA  : num  791 1059 743 1360 653 ...
##  $ E_MIN    : num  2678111 2673826 2675745 2680394 2674866 ...
##  $ E_MAX    : num  2681155 2678584 2679642 2686463 2678972 ...
##  $ N_MIN    : num  1234561 1235208 1238997 1230185 1238053 ...
##  $ N_MAX    : num  1238544 1239338 1243159 1236411 1241034 ...
##  $ E_CNTR   : num  2679300 2676800 2677800 2682900 2676400 ...
##  $ N_CNTR   : num  1235700 1236800 1241000 1233100 1239000 ...
##  $ Z_MIN    : num  533 440 502 524 475 463 412 386 432 383 ...
##  $ Z_MAX    : num  887 749 715 915 745 621 504 476 658 520 ...
##  $ Z_AVG    : num  685 529 581 673 560 557 454 414 506 440 ...
##  $ Z_MED    : num  673 502 583 653 543 559 454 410 493 439 ...
##  $ Z_CNTR   : num  700 490 544 610 503 573 431 400 471 453 ...
##  $ geometry :sfc_MULTIPOLYGON of length 2148; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:24, 1:2] 2678219 2678560 2678206 2678196 2678259 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:19] "GMDHISTID" "GMDNR" "GMDNAME" "BZHISTID" ...

Quick sf plots

# Quickly plot the geodata by selecting only the geometry
muni %>% select() %>% plot

# This will color the map by the features
# plot(muni)

# ggpplot2 way to visualize geo data
ggplot(muni) +
  geom_sf() 

Basic choropleth of Swiss municipalities by their canton

ggplot(muni) +
  geom_sf(aes(fill = KTNR))

Improved choropleth Swiss of municipalities by their canton

ggplot(muni) +
  # Hide the boundaries' line, discretize cantons
  geom_sf(aes(fill = as.factor(KTNR)), colour = NA) +
  # Hide colour scale
  scale_fill_viridis_d(guide = 'none') +
  # Remove graticule
  theme_void()

Get the votation results by municipality

muni_2022_05_15_results <- swissdd::get_nationalvotes(
  geolevel = "municipality", from_date = "2022-05-15", 
  to_date =  "2022-05-15", language = "FR")

# str(muni_2022_05_15_results)
#  muni_2022_05_15_results$name %>% unique()

ln_df <- muni_2022_05_15_results %>% 
  filter(id == 6550) %>% 
  select(canton_name:mun_name, jaStimmenInProzent, gueltigeStimmen)

Bind to geodata and plot

ln_df <- left_join(
  ln_df %>% mutate(mun_id = as.numeric(mun_id)), 
  muni %>% select(GMDNR:GRNR), 
  by = c("mun_id" = "GMDNR")
  # To ensure the join returns a sf object
) %>% st_as_sf()

ln_map <- ggplot(ln_df) +
  geom_sf(aes(fill = jaStimmenInProzent), colour = NA) +
  scale_fill_viridis_c(option = "E", limits = c(0, 100)) +
  theme_void() +
  theme(legend.position = "top") +
  labs(subtitle = "% de oui à la votation fédérale «Lex Netflix»")
ln_map 

Polish the map with cantons

# Recreate the cantonal boundaries with a spatial join

cant <- muni %>% 
  group_by(KTNR) %>% 
  summarise(geometry = st_union(geometry)) %>% 
  ungroup() 

plot(cant)

Better tooltip

# create more informative HTML tooltip
ln_df <- ln_df %>% 
  mutate(
    tp = str_c("<h3>", mun_name, "<h3>",
               "<i>", canton_name, "</i><br>",
               "<b>",
               prettyNum(jaStimmenInProzent, 
                         decimal.mark = ",", digits = 3),
               "% oui</b>"
               ) %>% 
      # ggiraph hack to ensure there's no single quote
      str_replace_all("'", "&#39;")
  )

ln_mapi <- ggplot(ln_df) +
  geom_sf_interactive(
    aes(fill = jaStimmenInProzent, tooltip = tp, data_id = mun_id), 
    color = NA) +
  scale_fill_viridis_c(option = "E", limits = c(0, 100)) +
  theme_void() +
  theme(legend.position = "top") +
  geom_sf(data = cant, size = 0.2, colour = "white", fill = NA)

sf cheatsheet